#include <xmmintrin.h>
#include <emmintrin.h>
#include <pmmintrin.h>
#include <tmmintrin.h>
#include <smmintrin.h>
#include <nmmintrin.h>
#include <immintrin.h>
#include <windows.h>
#include <iostream>
#include <stdlib.h>
#include <time.h>

using namespace std;

int main()
{
    long long freq,head,tail;
    srand(time(0));
    QueryPerformanceFrequency((LARGE_INTEGER *)&freq);
    int row = 0;
    for(int i = 1; i < 20;++i)
    {
        row += 100;
        QueryPerformanceCounter((LARGE_INTEGER *)&head);
        for(int q = 0;q<5;q++)
        {
            float ** mat;
            __m128 t1,t2,t3,t4;
            mat = new float*[row];
            for (int i = 0; i < row; i++)
            {
                mat[i] = new float[row];
                for (int j = 0; j < row; j++)
                    mat[i][j] = rand();
            }
            for (int k = 0; k < row; k++)
            {
                float temp[4] = { mat[k][k], mat[k][k], mat[k][k], mat[k][k] };
                t1 = _mm_loadu_ps(temp);
                float mid = mat[k][k];
                int suf = row - 1;
                for (int j = row - 1; j % 4 != 3; j--)
                {
                    mat[k][j] = mat[k][j] / mid;
                    suf = j;
                }
                for (int j = suf - 3; j % 4 == 0 && j >= k + 1; j -= 4)
                {
                    t2 = _mm_loadu_ps(mat[k] + j);
                    t3 = _mm_div_ps(t2, t1);
                    _mm_storeu_ps(mat[k] + j, t3);
                    suf = j;
                }
                for(int j = k + 1; j < suf; j++)
                    mat[k][j] = mat[k][j] / mid;
                mat[k][k] = 1.0;
                for (int i = k + 1; i < row; i++)
                {
                    float temp1[4] = { mat[i][k], mat[i][k], mat[i][k], mat[i][k] };
                    t1 = _mm_loadu_ps(temp1);

                    int end1 = row - 1;
                    for (int j = row - 1; j % 4 != 3; j--)
                    {
                        mat[i][j] -= mat[i][k] * mat[k][j];
                        end1 = j;
                    }
                    for (int j = end1 - 3; j % 4 == 0 && j >= k + 1; j -= 4)
                    {
                        t2 = _mm_loadu_ps(mat[i] + j);
                        t3 = _mm_loadu_ps(mat[k] + j);
                        t4 = _mm_sub_ps(t2, _mm_mul_ps(t1, t3));
                        _mm_storeu_ps(mat[i] + j, t4);
                        end1 = j;
                    }

                    for (int j = k + 1; j < end1; j++)
                    {
                        mat[i][j] -= mat[i][k] * mat[k][j];
                    }

                    mat[i][k] = 0;
                }
            }
            for (int i = 0; i < row; i++)
            {
                delete[] mat[i];
            }
            delete[] mat;
        }
        QueryPerformanceCounter((LARGE_INTEGER *)&tail );
        cout <<(tail-head)*200.0/freq <<endl;
    }
}

